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Abstract 

The research work outlined in the present note highlights the essential role played 
by the simulation procedures implemented by us on CINECA supercomputers to com- 
plement the mathematical investigations carried within our group over the past several 
years. The ultimate target of our research is the understanding of certain crucial fea- 
tures of the information processing and transmission by single neurons embedded in 
complex networks. More specifically, here we provide a bird's eye look of some ana- 
lytical, numerical and simulation results on the asymptotic behavior of first passage 
time densities for Gaussian processes, both of a Markov and of a non-Markov type. 
Several figures indicate significant similarities or diversities between computational 
and simulated results. 



1 Introduction to FPT models for the neuronal fir- 
ing 

The research work outlined here takes place within the framework of applied 
probability. Our aim is to describe the dynamics of the neuronal firing by 
modeling it via a stochastic process representing the change in the neuron 
membrane potential between each pair of consecutive spikes (cf., for instance, 
[2U\). In our approach, the threshold voltage is viewed as a deterministic 
function, and the instant when the membrane potential reaches it (i.e. when 
a spike occurs) as a first passage time (FPT) random variable. We shall 
focus our attention on neuronal models rooted on Gaussian processes, partially 
motivated by the generally accepted hypothesis that in numerous instances 
the neuronal firing is caused by the superposition of a very large number 
of synaptic input pulses which is suggestive of the generation of Gaussian 
distributions by virtue of some sort of central limit theorems. 

Let us first consider a one-dimensional non-singular Gaussian stochas- 
tic process {X{t),t >to} and a boundary S{t) G C^[to,+oo). We assume 
P{X{to) = xq} = 1, with xq < S'(to)) i-e. we focus our attention on the subset 
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of sample paths of X{t) that originate at a preassigned state xq at the initial 
time Iq. Then, 



r,.„ = inf {t : X{t) > Sit)}, xo < S{to) 

t>ta 



is the FPT of X{t) through S{t), and 

g{t\xo,to) = -^^P{T,,<t) (1) 

is its probability density function (pdf). 

Henceforth, the FPT pdf g{t\xo,tQ) will be identified with the firing pdf 
of a neuron whose membrane potential is modeled by X{t) and whose firing 
threshold is S{t). 

In order to include more physiologically significant features - such as a 
finite decay constant of the membrane potential, the presence of reversal po- 
tentials, time-dependent firing thresholds - and to refer to wider classes of 
inputs as responsible for the observed sequences of output signals released by 
the neuron, we also define the FPT upcrossing model. This is viewed as an 
FPT problem to a threshold, or boundary, S{t) for the subset of sample paths 
of the one-dimensional non-singular Gaussian process {X{t),t > to} originat- 
ing at a state Xq. Such initial state, in turn, is viewed as a random variable 
with pdf 

^^"^"^ xo < S{to) - e 



ls{xo,to)^l P{X{to)<S{to)-e} 

i 0, xo > S{to) - e. 

Here, e > is a fixed real number and /(xq) denotes the Gaussian pdf of 
X{to). Then, 

r|,)= inf{t:X(t)>5(t)}, 
is the e-upcrossing FPT of X{t) through S{t) and the related pdf is given by 

I \ d (^\ rS{to)-e 

g^uKAto) = olPi^Xo <i)= I g{t\xo,to)^eixo,to) dxo (t > to), 



oo 



where g{t\xo,to) is defined in (P). Without loss of generality, we set to = 

and Xo = and for this case we write g{t) = g{t\0, 0) and gu{t) = gu\t\0), for 
fixed values of e. 

The selection of one of the various methods available to compute the firing 
pdf's g{t) and gu{t) depends on the assumptions made on X(t). For diffusion 
processes (cf. 1 , 2 ) and for Gauss-Markov processes (cf. [S]) we have proved 
that the firing pdf is solution of a second kind integral Volterra equation. For 
generally regular thresholds we have designed, and successfully implemented. 
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a fast and accurate numerical procedure for solving such integral equation, and 
compared our approximations with those stemming out of standard numerical 
methods. Furthermore, in [8|, by adopting a symmetry-based approach, we 
have determined the exact firing pdf for thresholds of a suitable analytical 
form. 

Mathematical models based on non-Markov stochastic processes better 
describe the correlated firing activity, even though their analytical treatment is 
more complicated and only rare and fragmentary results appear to be available 
in literature. 

By using a variant of the method proposed in ^2]) for a zero-mean non- 
singular stationary Gaussian process differentiable in the mean square sense, 
a cumbersome series expansion for the conditioned FPT pdf density and for 
the upcrossing FPT pdf density has been obtained (Uni, [II])- In both cases 
we have succeeded to obtain numerically a reliable evaluation of the first term 
of the series expans ion. By comparisons of these results with the simulated 
firing densities ([£], UOl)) we have been led to conclude that this first term is 
a good approximation of firing pdf 's only for small times. 

The simulation procedures provide valuable alternative investigation tools 
especially if they can be implemented on parallel computers, (see f^). We 
point out that our simulation originates from Franklin's algorithm 15 . We 
have implemented it in both vector and parallel modalities after suitably mod- 
ifying it for our computational needs, for instance to obtain reliable approxi- 
mations of upcrossing densities (cf. [S], [S])- Thus doing, reliable histograms of 
FPT densities of stationary Gaussian processes with rational spectral densities 
can be obtained in the presence of various types of boundaries. 

For a thorough description of the simulation algorithm we refer the reader 

to m, [s]. 

We wish to stress that our endeavors strive to improve the simulation 
techniques, particularly within the context of stationary Gaussian processes 
for which alternative simulation procedures are being implemented and tested. 
Besides the use of algorithms based on special properties of the spectral den- 
sity, such as its being of a rational type, methods based on the spectral repre- 
sentation of the process and circulant embedding methods are presently under 
investigation. The aim is to design efficient algorithms to simulate Gaussian 
processes of a more general type. It must also be pointed out that within our 
research project the role of the simulation procedure is threefold: 

(i) to provide an investigation tool to explore the behavior of the firing 
density in a variety of different conditions; 

(ii) to permit us to evaluate reliability and precision of the results obtained 
via numerical and analytic approximations; 

(iii) to represent the only possible alternative whenever analytic and compu- 
tational methods are failing. 
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Figure 1: Plots refer to FPT pdf g(t) with /3 = 0.5 for a zero- mean Gaussian 
process characterized by correlation function jSjl in the presence of the boundary Q 
with d — 0.25. In Figure 1(a), the function g{t) has been plotted. The estimated 
FPT pdf g{t) with a = 10"^" is shown in Figure 1(b), with a = 0.25 in Figure 1(c) 
and with a = 0.5 in Figure 1(d). 



2 Markov and non-Markov models: an analysis by 
simulations 

In order to ana lyze how the lack of memory affects the shape of the FPT 
densities, in we have compared the behavior of such densities for Gauss- 
Markov processes versus Gauss non-Markov processes. 

Let us consider a zero-mean stationary Gaussian process X(t) with corre- 
lation function 



lit) 



-/3|t| 



cos (a t), a,pe]R'^ 



(2) 



whi ch is the simplest type of correlation of concrete interest for applications 
|22j . Furthermore, let us assume that the boundary is of the following type: 



Sit)=de-I^Ul 



2d2 



In 



11/ / 4d^ 

- + -Wl + 8exp(^--^^ 



(3) 



with d > 0. Due to the assumed correlation @, X{t) is not mean square differ- 
entiable (see for the details). Thus, the afore-mentioned series expansion 
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(see does not hold. However, specific assumptions on the parameter a 
help us to characterize the shape of the FPT pdf. 

We start assuming a = 0, so that the correlation function factorizes as 

^(t) = e-^^e-'^(*-^) (3 £ M+,0 <T <t. 

In such case X{t) becomes Gauss-Markov. As proved in [HI, for the Gauss- 
Markov process with covariance the FPT pdi g{t) in the presence of 
boundaries of type Q can be evaluated in closed form. Alternatively, g{t) can 
be numerically obtained by solving a non-singular Volterra integral equation 
(see |H]). In Figure 1(a) such density is plotted for /? = 0.5 and d = 0.25. 

Setting a 7^ in the Gaussian process X{t) is no longer Markov and 
its spectral density is given by 

^ ' Cj4 + 2 u;2 (^2 _ „2) + (^2 + a2)2 ' ^'^^ 

thus being of a rational type. Since in Q) the degree of the numerator is 
less than the degree of the denominator, it is possible to apply the simulation 
algorithm described in ^ in order to estimate the FPT pdf ^(t) of the process. 

The simulation procedure has been implemented by a parallel FORTRAN 
90 code on a 128-processor IBM SP4 supercomputer, based on MPI lan- 
guage for parallel processing, made available to us by CINECA. The number 
of simulated sample paths has been set equal to lO"^. The estimated FPT 
pdf's g{t) through the specified boundary are plotted in Figures l(b)-^l(d) 
for Gaussian processes with correlation function Q in which we have taken 

a = 10~^°, 0.25, 0.5, respectively. Note that as a increases, the shape of the 
FPT pdf g{t) becomes progressively flatter, while its mode increases. Further- 
more, as Figures l(a)-l(b) indicate, g{t) is very close to g{t) for small values 
of a. 



3 Asymptotic results 

The asymptotic behavior of the FPT densities for Gaussian processes as bound- 
aries or time grow larger has been studied in 9 , ,10 and 11 . Our analysis is a 
natural extens ion of some investigations performed for the Ornstein-Uhlenbeck 
(OU) process jlH] and successively extended to the class of one-dimensional 
diffusion processes admitting steady state densities in the presence of sin- 
gle asymptotically constant boundaries or of single asymptotically periodic 
boundaries (see ^4 and 17J. There, computational as well as analytical 
results have indicated that the conditioned FPT pdf is susceptible of an excel- 
lent non-homogeneous exponential approximation for large boundaries, even 
though these boundaries need not be very distant from the initial state of the 
process. To this aim, we have estimated such a density by generating the sam- 
ple paths of the Gaussian process through the parallel simulation algorithm 
implemented on the super-computer CRAY T3E of CINECA in order to over- 
come the outrageous complexity offered by the numerical evaluation of the 
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Figure 2: For a zero-mean stationary Gaussian process with covariance (01 and 
a = /3 = 1 in the presence of the periodic boundary S{t) =2 + 0.1 sin(27ri/3), the 
simulated FPT density g{t) is compared with the exponential density Ae""***, with 
A = 0.030386 in (a). The function Z{t) = g{t) exp{At} is plotted in (b). The 
asymptotic exponential approximation is plotted together with the simulated FPT 
density g{t) in (c). The asymptotic exponential approximation is plotted together 
with the simulated upcrossing FPT density gu{t) in (d). 



involved partial sums of the conditioned FPT pdf series expansion. Specifi- 
cally, we have considered the class of zero-mean stationary Gaussian processes 
characterized by damped oscillatory covariances (see 



7(i) 



'l^\t\ 



cos (a t) + sin(Q;|t|) 



(5) 



where a and /3 are positive real numbers. 

The results of our computations have shown that for certain periodic 
boundaries of the form 



S{t) = So + Asm{2TTt/Q), So,A,Q>0, 



(6) 



not very distant from the initial value of the process, the FPT pdf soon exhibits 
damped oscillations having the same period of the boundary. Furthermore, 
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starting from rather small times, the estimated FPT densities g{t) appears to 
be representable in the form 

m ^ Zit) e-^\ (7) 

where A is a constant that can be estimated by the least squares methods, and 
where Z{t) is a periodic boundary of period T (see, for example. Figure l2{a) 
and Figure Efb)). The goodness of the exponential approximation increases 
as the boundary is progressively moved farther apart from the starting point 
of the process. The more the periodic boundary is far from the starting point 
of the process, the more the exponential approximation improves. 

In ITP we have shown by rigorous mathematical arguments that as bound- 
ary Q moves away from the initial state of the process, the FPT pdf ap- 
proaches a non-homogeneous exponential density of the type 

~/i(t)exp|-^*/i(r)dT|, (8) 

where 



Mt) = ^^exp|-^j 



exp 



2(a2 + /32) 



" ./^(t)Erfc^ '^'^ 



2(a2 + /32)^'w (a'+p-^) 



(9) 



and p{t) = Asm{27rt/Q). (See, for example, Figure I21^c)). In 10 a similar 
result is proved for the upcrossing FPT density. (See Figure ISfd)). 

It should be stressed that the analytic and the simulation results consti- 
tutes only a preliminary step towards the construction of neuronal models 
based on non-Markov processes. Nevertheless, the unveiling of properties of 
the asymptotic behavior of FPT may turn out to be useful also for the descrip- 
tion of neuronal activities at small times whenever the intrinsic time scale of 
the microscopic events involved during the neuron's evolution is much smaller 
than the macroscopic observation time scale, or when the asymptotic regime is 
exhibited also in the case of firing thresholds not too distant from the resting 
potential, similarly to what was already pointed out by us in connection with 
the OU neuronal model 



4 An alternative approach 

Within the context of single neuron's activity modeling a completely different, 
apparently not well known, approach was proposed by Kostyukov et al. ( 16 ) 
in which a non-Markov process of a Gaussian type is assumed to describe the 
time course of the neural membrane potential. The model due to Kostyukov 
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Figure 3: For different clioices of in the interval [0.008,1.024], with threshold 
S{t) = — — t + 5, plot of the simulated is shown in (a) and plot of gu{t) for 
the OU-model in (c). Same in (b) and (d) for values of ^ in [2.048, 200]. 

(K-model) makes use of the notion of correlation time. Namely, let X{t) be 
a stationary Gaussia n p rocess with zero mean, unit variance and correlation 
function R{t). Then [21], 

= / \R(t)\ dT < +00 

Jo 

is defined as the correlation time of the process X{t). Under some assumptions 
on the threshold and by using a sort of diffusion approximation, Kostyukov 
works out a numerical evaluation q(t) to the upcrossing FPT pdf. This ap- 
proximation is obtained as solution of an integral equation that can be solved 
by routine methods. The relevant feature of this approach is that the unique 
parameter i? characterizes the considered class of stationary standard Gaussian 
processes. 

In 2 and in we analyzed this method pinpointing similarities and dif- 
ferences with respect to our models. Recently we have again made use of 
the K-model and compared the obtained results with those worked out by nu- 
merically solving the integral equation holding for Gauss-Markov processes in 
the case of the OU and Wiener models, see ,5^ for details, and with the results 
obtained via the simulations of Gaussian processes. A variety of thresholds 
and of -d values has been considered. For some values of between 0.008 and 
200, our results have been plotted in Figures 01 and |^ 
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Figure 4: For different cfioices of i? in tiie interval [0.008,1.024], with the same 
boundary as in Figure 01 plot of gu{t) for the Wiener model is shown in (a) and 
plot of q{t) for the Kostyukov-model in (c). Same in (b) and (d) for values of "d in 
[2.048,200]. 

Our investigations in this direction suggest that the validity of approxi- 
mations of the firing densities in the presence of memory effects by the FPT 
densities of Markov type is clearly depending on the magnitude of the corre- 
lation time. Hence, an object of our present research is the investigation of 
those models whose asymptotic behavior becomes increasingly similar as the 
correlation time 'd grows larger. 
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